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A Gutzwiller-type variational wave function is proposed for the neutron resonance mode in the 
cuprate superconductors. An efficient re- weighting technique is devised to perform variational Monte 
Carlo simulation on the proposed wave function which is composed of linearly superposed Gutzwiller 
projected Slater determinants. The calculation, which involves no free parameter, predicts quali- 
tatively correct behavior for both the energy and the spectral weight of the resonance mode as 
functions of doping. 
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INTRODUCTION 

The (tt, tt) resonance mode observed by neutron scat- 
tering in the cuparte superconductors is among the most 
prominent phenomena in these systemsJj,L^, ^S,, JJ. Below 
the superconducting transition temperature Tc, a sharp 
peak is observed in the spin fluctuation spectrum around 
(tt, tt). The mode is found to have close correlation with 
superconductivity of the system. For example, the mode 
energy, which is temperature independent, is found to 
scale linearly with Tc as a function of doping. At the 
same time, the intensity of the mode is found to have the 
similar temperature dependence as the superfluid den- 
sity. 

Much theoretical efforts have been devoted to the un- 
derstanding of the origin of the neutron resonance mode 
and its correlation with superconductivity^, |2, 0, S S 
llol |. Among these theories, the RPA-like theory, which 
takes the resonance mode as a spin-one bound state in the 
particle-hole channel(spin exciton) induced by the resid- 
ual attractive interaction between the Bogliubov quasi- 
particles in the superconducting state, is the most pop- 
ular. In this theory, the dynamical spin susceptibility is 
given by 



xo(q,t^) 



1 - c^(q)xo(q,^) ' 



(1) 



here xo(q7'^) is the bare spin susceptibility of the BCS 
superconducting state determined by both the band dis- 
persion and the gap function. U (q) is the phcnomenolog- 
ical RPA correction factor chosen to fit the experimen- 
tal data. In the RPA theory, the superconducting gap 
sets a natural energy scale for the resonance mode below 
which the mode is stable. At the same time, when the 
system approaches the antiferromagnetic ordering insta- 
bility under RPA correction, the mode will evolve into 
the Goldstone mode of the ordered state and its energy 
will approach zero. Thus, the energy of the resonance 
mode in the RPA theory is determined by both the mag- 
nitude of the superconducting gap and the strength of 
the antiferromagnetic correlation. 



Although the RPA theory can account for some as- 
pects of the neutron resonance mode, it is at the best, 
a phenomenological theory. The band dispersion of the 
quasiparticle and the RPA correction factor t/(q), on 
which the result of RPA calculation depends sensitively 
on, are subjected to fine tuning. An effort to combine 
the more microscopic slave Boson mean field theory and 
the RPA theory in the t — J model Q has resulted in too 
large a doping range(0 < x < 0.2, x is doping concen- 
tration) in which the system is unstable with respect to 
antiferromagnetic ordering. At the same time, neither 
the phenomenological RPA theory nor the RPA correc- 
tion on the slave Boson mean field theory of the t — J 
model respects the local spin sum rule of the t — J model, 
J dqdu)S{q,Ld) = (1 — a;)|, as a result of their neglect of 
the no double occupancy constraint of the t — J model. 



In this paper, we propose a variational description for 
the neutron resonance mode with a Gutzwiller projected 
wave function. Our approach can be taken as the gener- 
alization of the usual RPA theory into the Hilbert space 
satisfying the no double occupancy constraint of the t — J 
model. The approach has the virtual that it is parameter 
free: the RPA correction is automatically done through 
the variational procedure. We also devise an efficient al- 
gorithm to do Monte Carlo simulation on the variational 
wave function we proposed. Numerical calculation shows 
that our variational description of the neutron resonance 
mode capture its basic characteristics very well. 



The paper is organized as follows. In the next section, 
we introduce the variational ground state and present 
the result of the single mode approximation on which 
the result of our variational calculation would compare 
to. We then introduce our variational wave function for 
the neutron resonance mode and the numerical technique 
to do Monte Carlo simulation on it. Then we present our 
numerical results and offer a discussion on the results. 
Finally, we conclude this paper with some further prob- 
lems to be addressed in the future. 
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THE SINGLE MODE APPROXIMATION 

We take the t — J model as the basic model to describe 
the physics of high-Tc superconductors 

H = Ht + Ht>+ Hj, 

Ht - 



* XI (4<tCj,<t + /i.C.), 
<i,j>.cj 

t' (4.Jj.<T + h.c.), 



<<i,j>>,a 



1 



(2) 



Here Cj^a is the constrained electron operator satisfying 

the constraint E^cJ,<tC^<t < 1- J2<^,j> and 
represent the sum over nearest neighboring(NN) and 
next nearest neighboring (NNN) sites. Here we take 
-0.25 to describe hole doped system. The exchange 



t 

term is fixed at j 



The no double occupancy constraint X^ct ^1 o-^^-'^ — ^ 
is crucial for the spin dynamics of the system. With 
this constraint, the electron behave like mobile s — ^ 
spin rather than usual free electron. More quantitatively, 
the spin structure factor in the t — J model satisfies the 
following local spin sum rule 



J dqdLjS{q,uj) = {l — x) — , 



(3) 



in any physical state, here x is hole density. When the 
constraint is relaxed, as is done in slave Boson mean field 
theory or phenomenological RPA theory, the spin fluctu- 
ation would be strongly suppressed and no such sum rule 
would apply. 

To satisfy the local spin sum rule of the t — J model, 
the variational ground state, on which to construct the 
variational excitations, must respect the no double occu- 
panc y co nstraint. The Gutzwiller projected d-wave BCS 
statejlll 12, 13, which satisfy the no double occu- 
pancy constraint and for long has been known as an ex- 
cellent variational description of the ground state of the 
system, is the most natural choice for this purpose. 

Thus our variational ground state is given by 



(4) 



in which P/v is the projection operator into the subspace 
with N electrons and Pq is the Gutzwiller projection 
operator into the subspace of no double occupancy, a{i — 
j) — jig»k(ri-rj ) jg ^j^g j.gg^^ space wave function of the 

Cooper pah with ^ = Here, = V^+^ , 

and Ak are given by 

^k = — 2(cos kx + cos kj,) — At'^ cos k^: cos kj, — fj,^ 
Ak = 2Ai,(coska; — coskj,), (5) 



in which A^ are variational parameters to be de- 

termined by the optimization of the ground state energy 
with respect to the t — J model. We note t'^,fiv, are 
just variational parameters, rather than real NNN hoping 
term, real chemical potential and real superconducting 
gap. 

Now we construct the spin excitation on the variational 
ground state. As a first approximation to (tt, tt) reso- 
nance mode, we adopt the single mode approximation of 
the form 



(6) 



in which Sg is the creation operator of the spin density 
excitation at Q = (tt, tt). As the variational ground state 
satisfy the no double occupancy constraint, the spin exci- 
tation spectrum is guaranteed to obey the local spin sum 
rule. 

The single mode approximation is a good approxima- 
tion when the spin fluctuation spectrum is dominated 
by the contribution form the resonance mode. The exci- 
tation energy in the single mode approximation can be 
calculated in the standard way 



(*^I*"q> 



(7) 



where Eg denotes the ground state energy. Assuming 
that iJ|*g) = Egl^g), we have 



^ _l {^g\[SQ,[H,S^mg) 
^ 2 {^g\SQS+\^g) 

Using the commutation relation 



(8) 



(9) 



in which c] ^ is the constrained electron creation operator 
at site i and S'^' is spin lifting operator at site j, the mode 
energy in the single mode approximation can be shown 
to be given by 



^{^g\Ht\^g)-l{^g\Hj\^g) 

t^o = 



(10) 



It is interesting to note that the Hf term of the Hamil- 
tonian does not contribute when Q = (7r,7r). 

The mode energy calculated from Eq.(lO) is shown in 
Fig. 1. Two things are to be noted here. First, Eg in- 
creases monotonically with doping, consistent with ob- 
servation in the underdoped regime. Second, Eq ap- 
proaches zero at half filling in the thermodynamic limit. 
The monotonic increase of Eq is due to the rapid de- 
crease of the spin structure factor as a function of dop- 
ing which overcompensates the increase of the absolute 
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value of the kinetic energy and the exchange energy. The 
second is in fact a reflection of the Goldstone theorem. 
At half filling, the spin structure factor calculated from 
I'^g) scales super linearly with the number of lattice site 

{Ml = {'^g\SQS+\^'g)l{^g\^g) oc nI) as a result 
of the long range correlation induced by the Gutzwiller 
projection[l3j. However, the kinetic and exchange en- 
ergy is by definition proportional to Ns- Thus Eq must 
approach zero at half filling as A^^ — > oo. This is to be 
compared with the mean field prediction that AIq oc Ns 
at all doping which imply that Eq remains finite even at 
half filling. Thus the Gutzwiller projection in our wave 
function plays a key role to recover the correct trend of 
the mode energy as a function of doping. 

The mode energy calculated by the single mode ap- 
proximation is in fact the center of gravity of the spin 
fluctuation spectrum at Q = (tt, tt) and thus overesti- 
mates the energy of the resonance mode which lies at the 
bottom of the spin fluctuation spectrum. In fact, the sin- 
gle mode approximation we adopted above has nothing 
to say about the very existence of the resonance mode. 
This is especially clear at large doping when the mode en- 
ergy predicted by Eq.(lO) stretches into the particle-hole 
continuum. Thus, although the single mode approxima- 
tion gives correctly the trend of the mode energy as a 
function of doping, it is too crude to give a quantitative 
answer on both the mode energy and mode weight. 

It should be noted that the energy calculated from 
Eq.(lO) with the variational ground state is only a ap- 
proximation to the single mode approximation energy(as 
we have made the assumption that = Eg\'i>g)). 

Thus, although Eq.(IO) is by definition positive definite, 
the true single mode approximation energy can be neg- 
ative when the system become unstable with respect to 
magnetic ordering at Q = (tt, tt). We will encounter this 
situation below when we calculate the single mode ap- 
proximation energy directly from Eq.(7). 

THE PROJECTED SPIN EXCITON WAVE 
FUNCTION 

In the RPA theory, the resonance mode is interpreted 
as a spin exciton: a spin-one particle-hole bound state 
below the superconducting gap induced by the residual 
interaction between the Bogliubov quasiparticles. The 
wave function for the spin exciton can be generally writ- 
ten as 

I^Q)=E<^k7ilT7S-kTl'^-BCS), (II) 

k 

in which 0k describe the relative motion of the two quasi- 
particles 7j^-|. and 7Q_k| within the bound state, jd— BCS) 
denotes the BCS mean field ground state and 7j^-|, denotes 
the creation operator for the Bogliubov quasiparticles. 
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FIG. 1: The resonance energy determined from the single 
mode approximation as a function of doping. The inset shows 
the structure factor at Q = (tt, tt) (divided by the number of 
lattice sites, As ) as a function of doping. 

As we have shown, the mean field state fails to sat- 
isfy the local spin sum rule. For this reason, we project 
the spin exciton wave function Eq.(II) into the subspace 
of no double occupancy to construct a variational wave 
function for the resonance mode, with as the vari- 
ational parameter to be determined by optimization of 
energy. 

=PArPG|*Q) =5I'^k|k,Q), (12) 

k 

in which 

|k,Q) -PArPG7i[^7^_i^^|d-BCS). (13) 

We note that the wave function of the single mode ap- 
proximation can also be cast into the form of Eq.(12), 
i.e., 

Ivf^) = FArPc^+lrf-BCS) 

= ^PArPGcI+Q^^Ckald-BCS) 

k 

= E'^klkiQ), (14) 

k 

in which = Uk+Q^^fc- Here we have used the fact that 
commute with both and Pq- 
The variational parameters 0k can be determined by 
minimizing the energy of \^q) with respect to the t — J 
Hamiltonian. The variational energy is given by 

^ {^Q\Ht_j\^Q) ^ Ek.k- '/'k-ffk.k'^k' 

Ek,k'</'kOk,k'0k'' ^ ^ 

in which 

i?k,k' = (k,Q|i/t-,/|k',Q) (16) 
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is the matrix element of the t 
basis |k, Q) and 



J Hamiltonian in the 



Ok,k' = (k,Q|k',Q) 



(17) 



denotes the overlap integral between these non- 
orthogonal basis function(note that 7k|7Q-kTl'^ ^ BCS) 
form a orthogonal basis set before the Gutzwiller projec- 
tion). 

The problem of minimizing Eq respect to the set of 
variational parameters 0k now reduces to solving the fol- 
lowing generalized eigenvalue problem 



-ffk,k'0k' = X! ^k,k'0k 



(18) 



It is easily seen that the optimized energy Eq is given by 
the lowest eigenvalue A of the above generalized eigen- 
value problem. 

The above optimization procedure can also be inter- 
preted as re-diagonalizing the t — J Hamiltonian in the 
subspace spanned by the set of non-orthogonal basis func- 
tion |k, Q). With this understanding in mind, we can 
even construct variationally the full spin fluctuation spec- 
trum as follows 



S{Q,Lo)^M^ 



Q X! 

n.k.k' 



in which denotes the n-th eigenvector of the gener- 
alized eigenvalue problem with the eigenvalue En, Eg 
denotes the variational ground state energy(here we as- 
sume that (f>k and 0° are so chosen that both |5'q) 
and I^Pq) are normalized). As cj)]^ forms a orthonor- 
mal basis with respect to the overlap matrix Ok.k', *-e-, 
Ek k' K*Ok,k'<t>]^ = Sn,m, we have 



fdcjS{Q,u;) = M^^|5]0rOk,k'0^f 

•' n k,k' 

= M^^0O*Ok,k'0^. =m2, (20) 

k,k' 

in which we have used the fact that 

ECOk..k0L'-4,.k„ (21) 

n,k 

which can be derived from the orthonormality of 
the eigenvectors of the generalized eigenvalue problem 
Eq.(18), 



E0rok,k'C = '5„ 



(22) 



k,k' 



Thus our variational construction of the spin fluctuation 
spectrum respects the local spin sum rule of the t — J 
model. 



Our variational scheme for the resonance mode has the 
advantage that it involves no tunable parameter: the pa- 
rameters in l^g) is determined by optimizing the ground 
state energy and the spin exciton wave function 0k is de- 
termined by solving the generalized eigenvalue problem. 
The RPA correction is thus automatically taken into ac- 
count in our formalism. 

From the above discussion, we know the single mode 
approximation can be taken as a special case of the 
projected spin exciton wave function(with 0k = 0° = 
UkWQ-k )■ Thus by construction, the resonance energy 
calculated from Eq.(18) should be lower than that cal- 
culated from single mode approximation. In fact, the 
mode energy calculated from the single mode approxi- 
mation gives the center of gravity of the spin fluctuation 
spectrum, while the resonance mode lies at the bottom 
of the spectrum. As the single mode approximation al- 
ready reproduces the correct trend for the mode energy 
as a function of doping in the underdoped regime, we can 
even expect our variational scheme to produce quantita- 
tively reasonable result. 



THE RE- WEIGHTING TECHNIQUE 

To calculate the energy of the resonance mode, we 
should first evaluate the matrix element of the Hamil- 
tonian in the set of strongly correlated basis func- 
tions |k, Q), -ffk.k' and the overlap matrix element Ok.k'- 
This can be done in principle by the variational Monte 



Ok,i 



Carlo(VMC) method. For example, to evaluate — 

^^''^{'k^Q\k^Q)^^ ' '^'^^^ expand |k, Q) in an orthogonal 
basis i.e., 



|k,Q) = EVk(i?^)|i?^), 



(23) 



where ip{Ri) is the wave function of |k, Q) in this basis. 
Then we have 



Ok,k 



in which 



(24) 



(25) 



Then we sample the basis space \Ri) with the weight 
li^kiRi)]'^ and do the sum with the standard VMC 
technique [13] ■ 

The above procedure, though straightforward, is very 
inefficient. In our problem, there are Hamiltonian 
matrix elements i?k,k' and iV^ overlap matrix elements 
Ok.k' to be evaluated. For lattice of reasonable size, say, 
14 X 14, the number of the matrix elements to be eval- 
uated would exceed 40000 even take into account the 
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Hermitian property of iJk.k' and Ok,k'- This is very time 
consuming. At the same time, the naive approach has 
the drawback that it involves large statistical error in 
the simulation. This can be seen as follows, 



O 



k,k' 



(26) 



Thus, when the node of ■0k(^i) and ip\^i{Ri) do not co- 
incide with each other, wc will run into trouble when 
we sample ip]^{Ri) around its node as the fluctuation of 
'^Ur'} becomes large. A way to reduce the statistical er- 
ror caused by the fluctuation of ^^'(^ff is to sample the 
combined weight W(iii) = |V'k(-Ri)|^ + |V'k'(-Ri)|^ rather 
than |Vk(i?i)P, 



k.k' 



o 



k,k' 



Ok.l 



Ok.k + Ok'.k' 



/ 



o 



k,k 



o 



k.k - 



Ok' 



(27) 



Now the calculation of 



Ok,k' 



can be done as 



Ok,k + Ok', k' 



(28) 



The combined weight samples symmetrically between 

tp\^(Ri) and V-'k'(-Ri) and avoids the fluctuation caused 
by their uncommon nodes. Thus the statistical error is 
much reduced. 

The above technique can be easily generalized to cal- 
culate all the overlap matrix elements Ok,k'- Here 
we sample the combined weight of all the Ns |'0k(-Ri)p: 
WiRi) = X^k \i'kiRi)\'^- The calculation is done as fol- 
lows: 



' W{Ri) 



EkOk.k 



Eft MRi) 



To sample yV{Ri), we note that 

w{Ri) = \Aoim^T. 



V'ko(i?.)' 



and 



^ I i>^{Ri) 12 



(29) 



(30) 



(31) 



in which ^ko is one basis function arbitrarily chosen from 

the Ns basis functions. From this transformation, wc sec 
all we need to calculate in order to evaluate ^ 



the Ns ratio between basis functions 



MR^) 



IS 



As the 



different basis functions are all Slater determinant differ- 
ing with each other by at most in a pair of quasiparticle 

excitations, such ratio is easy to calculate using the in- 
verse updating technique for Fermion determinant. More 



importantly, the calculation of all the A'^^ overlap matrix 
elements can be done in a single Monte Carlo simulation: 
the algorithm is highly parallelized. 

The calculation of the Hamiltonian matrix elements 
can be done similarly. We have 



i?ki,k2 ^Ri ^^y^^i w(R,) 



EkOk,k 



(32) 



where Hi,],,{Ri) = Y.R,{Ri\Ht-j\Rj)i>kARj)- Follow- 
ing the same reasoning, we arrive at 



i*l-n\ TT I rr>\ / V'ki(fli) x* Hjik^jRi) 



W{Ri) 



(33) 



Thus the calculation of the Hamiltonian matrix elements 
involves the evaluation of the ratio _ Ti^g calcu- 

lation of this ratio, though numerically more demanding, 
is still highly parallelized. 

Thus the re-weighting technique developed here not 
only reduce considerably the statistical error involved in 
the Monte Carlo simulation, but also highly parallelize 
the calculation of the overlap and Hamiltonian matrix 
elements, reducing their calculation from the order of N^ 
to a single Monte Carlo simulation. 



NUMERICAL RESULTS 

To calculate the energy of the neutron resonance mode, 
we first optimize the variational parameters i^, iiy, A„ as 
a function of doping for the ground state. The calculation 
is done on a 14 x 14 lattice with pcriodic-antipcriodic 
boundary condition. Wc choose j = and j — —0.25 
in the t — J model to describe a hole doped cuprate. 
The results of the optimized variational parameters as a 
function of doping arc shown in Fig. 2 and Fig. 3. 

We note the superconducting region determined by the 
variational approach for the t — J model is considerably 
larger than that observed in experiments. In Fig. 3, we 
plot the off diagonal long range order(ODLRO) calcu- 
lated from the optimized variational ground state. The 
ODLRO is defined as 



(34) 



in which = (c,+j;,|Ci_| - Q+2;jc,_x) - (cj+y,|Cj_| - 
Ci+y,|Ci,l) is the Cooper pair annihilation operator at site 
i, Rm is the largest distance on the finite lattice. We find 
the ODLRO reaches its maximal around x w 26%. Wc 
will take this doping concentration as an estimate of the 
location of the optimal doping in the following discussion. 

After obtaining the variational parameters of the 
ground state we are now ready to calculate the Hamil- 



FIG. 2: The optimized variational parameters as a function 
of doping. 
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FIG. 3: The optimized pairing parameters as a function of 
doping and the off diagonal long range order calculated from 
the variational ground state. 



tonian matrix elements -ffk.k' and the overlap matrix el- 
ements Ok,k'- This is the most heavy part of our nu- 
merical calculation. In our calculation, we have sampled 
more than two and a half million configurations with the 
weight W(i?i). The accept ratio is tuned to be i. The 
statistical error is found to be smaller than the fluctua- 
tion caused by the finite size effect in our calculation. 

When we get the matrix elements i?k,k' and Ok,k', we 
solve the generalized eigenvalue problem Eq.(18) and cal- 
culate the physical quantities that we are interested in. 
First we show the spin fluctuation spectrum calculated 
in this variational approach in Fig. 4. We see the spec- 
trum consists of both a coherent peak and a continuum 
of incoherent spin fluctuation. The coherent peak at the 
bottom of the spectrum is nothing but the neutron reso- 
nance mode in our variational description. 



FIG. 4: The spin fluctuation spectrum at g = (tt, tt) as a 
function of excitation energy and doping determined from the 
variational calculation. The calculation is done on a 14 x 14 
lattice. The delta function peaks of the spectrum are broad- 
ened into Lorentzian peaks with a width of O.lt. 
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FIG. 5: The mode energy determined by our variational ap- 
proach as a function of doping as compared with the result of 
single mode approximation(SMA). 



Fig. 5 shows the mode energy as a function of doping. 
We find the mode energy becomes negative below a crit- 
ical doping around x — 7.5%. A negative excitation en- 
ergy indicates magnetic instability of the system. The 
critical doping so determined is close but sightly lower 
than that determined by assuming directly a magnetic 
order in the variational ground state, which is about 
10% IJ]. In whole doping range in which the mode has 
nonzero spectral weight, the mode energy is a monotoni- 
cally increasing function of doping and reaches about 0.3i 
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before it loses its weight and merges into the particle-hole 
continuum at about x = 29%. If we take t = 0.25eV as 
is usually done in the literature, we get the maximum 
of the mode energy to be about 75meV, about a fac- 
tor of 1.8 larger than that observed in optimally doped 

YBa2Cu306.93^. 



come from the long range correlation of spin at half filling 
and indicates that the resonance mode can be connected 
smoothly to the Goldstone mode in the magnetic ordered 
state. 



DISCUSSION 
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FIG. 6: The relative spectral weight of the resonance mode 
as a function of doping. The inset shows the result of the 
absolute spectral weight. 

In our theory, the mode energy increases monotoni- 
cally with doping. This is consistent with experiments 
in the underdoped regime but may have inconsistency 
with experiment in the overdoped regime, where experi- 
ment reported evidence of a weak resonance mode with 
an energy sightly lower than that of the optimally doped 
system [4j. If we take the doping at which the ODLRO 
reaches maximum(x = 26%) as the optimal doping, then 
the resonance mode will survive in the slightly overdoped 
regime with an energy higher than that of the optimal 
doped system. The would imply a breakdown of the 
linear scaling between the mode energy and in the 
slightly overdoped regime. 

Fig. 6 shows the relative and absolute spectral weight 
of the resonance mode calculated from our theory. The 
relative spectral weight is defined as the proportion of the 
mode intensity to the total spectral weight at Q = (tt, tt). 



W = 



I(*°qI*Q>P 



(35) 



The absolute spectral weight is defined as the product of 
the relative spectral weight and the spin structure factor 
divided by Ng. The absolute spectral weight decreases 
rapidly with doping as a result of the decrease of both 
the magnetic structure factor and the relative spectral 
weight. The relative spectral weight of the resonance 
mode decrease from unity at half filling down to zero at 
about X = 29%. The unity of the relative spectral weight 



In this paper, we proposed a variational theory for the 
neutron resonance mode in the cuprate superconductors. 
Our theory has the virtue that it involves no free pa- 
rameter and thus has much larger predictive power than 
the phenomenological RPA treatment of the spin fluctu- 
ation. In our theory, the RPA correction is automatically 
taken into account by the variational procedure which re- 
duces to solving a generalized eigenvalue problem. More 
importantly, our variational approach builds in the no 
double occupancy constraint and thus satisfies the local 
spin sum rule of the t — J model. This is of vital impor- 
tance for a correct description of the spin dynamics of 
cuprates. Our approach has the further advantage that 
it provides a physical transparent understanding of the 
resonance mode as a spin exciton in the physical subspace 
of no double occupancy. 

Our approach can also be taken as diagonalizing the 
Hamiltonian in a truncated subspace with the same quan- 
tum number as the excitation discussed and can be used 
to calculate the full spectrum rather than the coherent 
excitation only. It is important to note this truncated 
subspace exhaust the spectral weight for the relevant sum 
rule, indicating the relevance of the spectrum calculated 
in this way. An effort to apply the current approach 
to calculate the spin fluctuation spectrum at momentum 
other than (7r,7r), namely the incommensurate spin fluc- 
tuation spectrum, is now under investigation. 

We have also devised a very efficient re- weighting tech- 
nique to tackle the numerical problem of simulating the 
variational wave function composed of N linearly super- 
posed Slater determinants. The key for the efficiency of 
the algorithm is the observation that the Slater determi- 
nants involved in our wave function differ with each other 
by at most a pair of quasiparticle excitations. Obviously 
this technique can be applied in a much larger literature 
than simulating the physics of cuprates. 

As we have found in the single mode approximation, 
the center of gravity of the spin fluctuation spectrum in- 
creases monotonically with doping as a result of the rapid 
decrease of the spin structure factor at Q = (7r,7r). We 
find the resonance mode, which lies at the bottom of the 
spin fluctuation spectrum, inherits this monotonic behav- 
ior, probably for the same reason. As we have mentioned 
above, this may have potential conflict with the report of 
weak resonance mode in the overdoped sample with an 
energy lower than that of the optimal doped systemjl] 
and would imply the breakdown of the Tc — Er linear 
scaling in the overdoped regime. 
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In our theory, the monotonic increase of the mode en- 
ergy eventually cut off at about twice the maximal su- 
perconducting gap when the mode transfers all of its 
weight into the particle-hole continuum. To exhibit a 
non-monotonic behavior before merging into the particle- 
hole continuum, it is necessary for the superconducting 
gap to decrease faster with doping than that predicted 
by the present variational calculation. This is not at all 
impossible. However, since there is no generally reason 
to believe the Tc — Er linear scaling to hold in the over- 
doped regime and the mode weight become very small in 
the overdoped regime, we think the mode energy in the 
overdoped regime is a problem subjected to fine-tuning. 

The mode energy at the optimal doping(a: = 26%) as 
calculated from our approach is about 0.3i and is a fac- 
tor of 1.8 larger than that observed in optimally doped 
YBa2Cu30Q,g3. At the same time, the variational the- 
ory predicts a considerably larger value of optimal doping 
than observed in experiments(a; = 16%). It is likely that 
these two problems and the problem of the Tc — Er lin- 
ear scaling to hold in the overdoped regime are related 
with each other. However, it is not clear to what extend 
should we attribute these disagreements with experiment 
solely to the limitation of the variational approach we 
have adopted rather than the intrinsic properties of the 
t — J model. 

Finally, we discuss the relation between the resonance 
mode and the superconductivity. As many other the- 
ories of the resonance mode[i, 0, i, our theory also 
predicts that the resonance mode becomes stronger and 
stronger with decreasing doping and evolves smoothly 
into the Goldstone mode of the ordered state at half fill- 
ing. On the other hand, in the phenomenological S0(5) 
theory of cuprates[loj|. in which the mode is understood 
as a pseudo-Goldstone mode accompanying breaking of 
an SO (5) symmetry between the d-wave superconduct- 
ing order and the antiferromagnetic order, the mode in- 
tensity is predicted to be proportional to the off diago- 
nal long range order(ODLRO) of the system. Thus the 
theory predicts that the mode intensity should decrease 
when we increase temperature or decrease doping and 
disappear out of the superconducting dome. On the ex- 
perimental side, the resonance mode is observed to loss 
weight with increasing temperature and to disappear in 
the normal state at optimal doping. In slightly under- 
doped sample, a broadened and weak signal is observed 
above T^. With further decrease of doping, the normal 
state signal becomes stronger and stronger and the en- 
hancement due to superconductivity becomes less and 
less prominent P, 0, S]- 

In our theory, we only consider the zero temperature 
case. At finite temperature, both quasiparticle exci- 
tations and collective fluctuation of the superconduct- 
ing order parameter will be thermally excited. The 
latter excitation is believed to be especially important 
around the superconducting transition point. We be- 



lieve both of these thermal excitations are responsible 
for the decrease of the mode intensity with increasing 
temperature. The doping dependence of the mode in- 
tensity is more subtle. Here the relevant question is 
why the superconductivity-related enhancement of the 
mode intensity becomes smaller and smaller with de- 
creasing doping. In the very underdoped regime, the 
electron correlation(Mott physics)is greatly enhanced. 
We believe the Mott physics is at work in reducing the 
superconductivity-related enhancement of the mode in- 
tensity in the very underdoped regime. It is interesting 
to note that in the RVB picture, the electron correla- 
tion effect manifests itself as quantum fluctuation of the 
superconducting order parameter. 

In all, apart from some subtle issues mentioned above, 
our variational approach capture the gross feature of the 
neutron resonance mode and provides the first truly mi- 
croscopic understanding of this important phenomena in 
the high-temperature superconductors. As a by product 
of this research, we developed a very efficient re- weighting 
technique to simulated wave function composed of N 
linearly superposed Slater determinants, an algorithm 
whose potential application is obviously far beyond the 
high-Tc issue. 
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